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ABSTRACT 


This  short  interim  report  summarizes  the  work  completed  during  the 
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Research  as  part  of  Grant  No.  AFOSR77-3336.  The  research  deals  with  the 
broad  topics  of  initiation,  combustion  and  transition  to  detonation  in 
homogeneous  and  heterogeneous  reactive  mixtures.  One  specific  area  deals 
with  analytical  and  experimental  work  directed  to  direct  initiation  of 
detonation  by  a nonideal  blast  wave  in  chemically  sensitized  reactive 
fuel-air  clouds.  The  other  specific  topic  involves  the  hydrodynamic  model- 
ling of  ignition  and  flamespreading  in  granular  energetic  solids  to  predict 
the  potential  for  def lagration-to-detonation  (DDT) . 
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INTRODUCTION 

The  direct  initiation  of  detonation  by  a localized  source  such  as  a 
laser  pulsed  spark,  capacitance  spark,  exploding  wire,  or  exploding  high  ex- 
plosive charge  is  well  understood  at  the  present  time.  Initiation  by  the 
use  of  localized  chemical  accelerators  is  much  less  understood.  Here,  there 
is  still  the  question  of  whether  one  can  generate  an  initiating  shock  wave 
by  using  the  proper  distribution  of  a chemical  accelerator  in  a source  re- 
gion. Also,  nonlinear,  two-dimensional  initiation  behavior  is  not  understood. 
The  proposed  work  will  encompass  an  experimental  program  in  which  nonlinear 
initiation  behavior  will  be  studied  in  systems  of  direct  interest  to  the  Air 
Force,  and  a theoretical  program  in  which  the  basic  mechanisms  of  initiation 
by  localized  accelerator  concentrations  will  be  explored. 

It  is  expected  that  these  developments  will  justify  the  use  of  direct 
testing  required  to  discover  the  initiation  behavior  of  a wide  variety  of 
fuel-accelerator  combinations.  It  therefore  has  direct  applicability  to 
Phase  III  FAE  development  as  well  as  the  studies  of  fuel  tank  vulnerability. 
Furthermore,  the  program  will  have  application  to  understanding  the  nature  of 
nonideal  explosions  that  occur  during  explosive  release  accidents. 

The  DDT  phenomenon  in  granulated  propellant  or  explosives  involves  a 
series  of  complex  transient  processes  that  are  not  well  understood  at  the 
present  time.  It  is  hypothesized  that  the  normal  burning  process  of  the  solid 
propellant  is  disturbed  by  an  abnormality  such  as  a crack  in  the  propellant 
grain.  This  abnormality  generates  regions  of  porous  propellant  which  can  be 
ignited  locally,  causing  a pressure  buildup  and  formation  of  a weak  shock. 

If  detonation  is  to  be  excited  following  this  ignition,  it  is  necessary  to 
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ensure  a sufficiently  rapid  pressure  buildup.  In  the  case  of  porous  pro- 
pellant, this  may  be  achieved  as  a result  of  the  penetration  of  gaseous  com- 
bustion products  into  the  interior  pores  of  the  solid,  which  leads  to  the 
disturbance  of  surface  burning  conditions.  Thus,  in  this  case,  heat  transfer 
by  conduction  is  replaced  by  convective  heat  transfer.  Subsequent  accelera- 
tion of  ignition  (flame)  fronts  begins  and  pressure  waves  are  generated  which 
become  shocks.  These  shocks  cause  large  local  over-pressurization  and  often 
change  into  a detonation. 

To  analyze  this  phenomenon,  the  reactive  two-phase  (solid,  gas)  con- 
servation equations  of  continuity,  momentum  and  energy  must  be  solved  along 
with  many  constitutive  relations  to  account  for  heat  transfer  interaction, 
pressure  losses  through  the  aggregate,  ignition  criteria,  unsteady  burning 
rates,  etc.  It  is  only  through  solutions  of  such  a fluid -mechanics  model  that 
one  can  develop  criteria  to  specify  the  conditions  under  which  the  burning 
propellant  is  susceptible  to  transition  to  detonation.  Considerable  effort 
is  involved  In  developing  numerical  integration  schemes  to  adequately  handle 
the  prediction  of  strong  shock  waves  in  such  transient  flows. 
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SUWIARY  OF  PROGRESS 

During  this  period,  study  continued  on  phenomena  associated  with  ini- 
tiation, combustion,  transition  to  detonation  (DDT) and  attenuation  in  homo- 
geneous and  heterogeneous  reactive  media.  A modified  Lagrangian  time-depen- 
dent finite  difference  (Oppenheim  (CLOUD))  program  was  used  to  study  direct 
initiation  of  detonation  by  a nonideal  blast  wave  in  chemically  sensitized  re- 
active fuel-air  clouds.  Modeling  of  the  Arrhenius  kinetics  throughout  the  ex- 
plosion region  in  such  a manner  as  to  also  satisfy  computer  stability  require- 
ments and  reasonably  short  run  times  has  been  completed.  Testing  of  different 
sensitizer  concentration  profiles  in  the  source  region  has  been  started.  The 
shock  tube  study  of  initiation  using  a 15°  ramp  on  a splitter  plate  at  the 
back  wall  has  been  completed. 
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Studies  of  flame  acceleration  processes  by  obstacles  placed  between 
two  flat  plates  has  been  started  as  have  studies  of  flame  acceleration  pro- 
cesses by  obstacles  in  a layered  combustible  mixture  which  is  open  at  the  top. 

The  fluid  dynamics  needed  to  predict  def lagration-to-detonation  tran- 
sition in  granulated  beds  of  high  energy  solid  propellant  is  still  being 
modeled.  We  have  determined  the  sensitivity  of  the  inter -phase  viscous  drag 
and  heat  transfer  on  the  predictions  of  the  flame spreading  rate  in  such  packed 
beds.  Some  of  the  results  of  this  work  were  published  in  the  17th  Symposium 
(International)  on  Combustion.  The  paper  is  included  as  Appendix  A.  Also 
completed  during  this  past  year  were  (a)  a detailed  design  and  fabrication  of 
an  experiment  that  will  provide  the  necessary  pressure  drop  relation  through 
packed  beds  of  solid  particles,  as  would  be  needed  under  the  high  pressure 
convective  flow  conditions,  and  (b)  new  numerical  integration  schemes  that 
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A SYNOPSIS  OF  WORK  RELATED  TO 
FUEL/AIR  INITIATION  AND  BLAST  WAVES 


PROGRESS 

In  the  theoretical  initiation  program  we  have  adapted  the  Oppenheim 
CLOUD  program,  which  is  a constant  time  step  Lagrangian  finite  difference 
scheme,  so  we  can  examine  the  question  of  what  type  of  distribution  of  ac- 
celerator must  be  present  in  a spherical  source  region  to  generate  a shock 
wave  in  the  surrounding  region  of  sufficient  strength  to  cause  direct  initia- 
tion of  detonation  in  the  surrounding  region.  For  our  now  completed  first 
case,  the  distribution  of  accelerator  had  a uniform  central  core  surrounded 
by  a decreasing  concentration  given  by 

F(R)  = F1{cosC3tt6)  - 9.0  cos(tt<$]/16.0 

where  F^  is  the  concentration  in  the  central  core, 

R - R 

6 = ^ 5—  for  the  range  R.  1 R 1 R , 

K ” K lO 

O 1 

R,  is  the  radius  of  the  central  core,  and  R is  maximum  radius  at  which 
1 o 

accelerator  is  present.  In  this  first  case  we  made  Rj  = 0.2  Rq. 

We  first  attempted  to  use  an  Arrhenius  kinetic  law  throughout  the  ex- 
plosion region,  but  we  found  the  use  of  this  law  led  to  rates  of  energy  ad- 
dition becoming  so  large  that  the  numerical  stability  criterion  required  very 
small  time  steps  and  consequently  the  computation  time  became  excessive.  We 
then  modified  the  law  by  imposing  a maximum  rate  which  could  not  be  exceeded 
and  in  addition  imposed  a minimum  allowable  time  step. 

The  program  now  runs  well  with  no  instabilities.  After  verifying  that 
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we  could  run  the  program  adequately,  we  have  modified  the  Arrhenius  rate  law 
such  that  at  high  temperature  it  adequately  models  the  methane  oxidation  pro- 
cess and  below  1000°K,  has  an  extremely  high  effective  activation  energy  such 
that  the  reaction  just  cannot  occur  at  all  at  room  temperature.  The  high  and 
low  temperature  regimes  are  flared  together  smoothly  so  that  there  is  no  dis- 
continuity in  a plot  of  log  k versus  1/T.. 

We  have  successfully  obtained  detonation  initiation  for  cases  where 
the  central  core  is  highly  reactive  out  to  0.2  of  its  full  radius  and  is  flared 
as  a cosine  function  of  reactivity  out  to  the  edge.  We  have  also  successfully 
obtained  detonation  initiation  in  the  surrounding  gas  for  the  case  when  the 
flaring  extended  over  only  six  cells  (from  cell  46  through  cell  51.)  However, 
if  we  use  a bursting  sphere  for  the  central  region  we  do  not  get  initiation  in 
the  surroundings.  It  appears  from  the  output  of  the  calculation  that  the 
reason  for  this  is  that  when  the  central  region  reacts  very  rapidly,  inertial 
containment  of  the  very  center  of  the  region  causes  it  to  explode  more  rapidly 
than  the  neighboring  region.  The  resulting  compression  wave  has  temperatures 
within  it  which  are  higher  when  it  reaches  the  core  edge  than  the  temperature 
is  just  behind  the  initial  shock  wave  created  by  the  bursting  sphere  and  thus 
initiation  can  occur.  Further,  the  minimum  required  energy  release  rate  is 
greater  for  the  short  reactivity  transition  zone  than  it  is  for  the  wider  tran- 
sition zone.  In  addition,  there  is  the  possibility  that  extremely  large  values 
of  core  energy  release  rates  may  cause  the  compression  wave  to  change  into  a 
fully-developed  shock  wave  before  the  wave  reaches  the  edge  of  the  core.  If 
the  shock  wave  is  not  strong  enough,  then  the  temperature  may  be  insufficient 
to  cause  initiation  outside  the  core.  The  implication  of  these  results  is  that 
sufficiently  rapid  reaction  of  a central  region  in  the  cloud  can  trigger 
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detonation  irrespective  of  its  shape.  This  is  an  important  consideration 
because  it  means  that  one  need  not  be  too  fussy  about  the  shape  of  the  core 
reactive  region  in  order  to  get  detonation  provided  the  energy  release  rate 
in  the  reactive  region  is  within  an  appropriate  range.  Shaping  does  affect 
the  minimum  required  energy  release  rate,  however,  with. a larger  Value  needed 
for  a short  transition  region  than  is  needed  for  a wide  transition  region  for 
the  same  initial  shape. 

The  nonlinear  initiation  experiments  using  a 15°  ramp  at  the  back  wall 
of  the  shock  tube  were  quite  unsuccessful.  For  some  reason  we  were  never  able 
to  get  reasonable  smoke  track  records  and  because  of  the  two-dimensionality  of 
the  flow  we  made  no  attempt  to  take  streak  records  of  the  phenomenon.  That 
program  has  been  terminated. 

A program  has  been  started  to  study  the  effect  of  obstacles  in  the  path 
of  the  flame  on  flame  acceleration  processes.  Initially  we  used  two  flat 
plates  placed  about  2 inches  apart  and  blew  very  large  soap  bubbles  of  methane- 
air  stoichiometric  mixtures.  These  were  ignited  centrally  as  cylindrical 
bubbles  and  in  some  preliminary  experiments  a variety  of  obstacles  were  placed 
in  their  path.  The  obstacles  were  cylindrical  obstacles  surrounding  the  ig- 
nition source.  Obvious  accelerations  occurred  because  of  the  increase  in 
noise  level  when  the  obstacles  were  in  place.  However,  instrumentation  was 
not  developed  to  record  any  of  the  acceleration  phenomena.  We  have  recently 
again  shifted  our  experimental  emphasis.  We  are  now  planning  to  run  an  ex- 
periment in  which  a two-dimensional  flame  is  allowed  to  pass  over  two-dimen- 
sional objects  between  two  glass  plates.  Schlieren  observations  will  be  used 
to  record  the  flame  acceleration  processes.  This  experiment  is  different  from 
the  previous  one  in  that  the  chamber  which  contains  the  combustible  mixture 


will  be  open  at  the  top  and  will  be  only  partially  filled  with  combustible 


mixture  by  using  a layering  process.  The  point  here  is  that  free  clouds  are 

open  at  the  top  and  are  heavily  layered.  Thus  we  feel  it  is  more  realistic 

' 

to  examine  acceleration  by  obstacles  under  conditions  which  more  closely 
match  those  in  real  clouds.  Various  rectangular  obstacles  will  be  placed 
across  the  chamber  at  different  spacings  and  at  different  heights  relative 
to  the  height  of  the  layered  portion  of  the  cloud.  Repetitive  photographs 
using  the  schlieren  system  will  be  taken  to  observe  both  the  motion  of  the 
cloud  and  the  motion  of  the  flame  inside  the  cloud.  Optimum  conditions  for 

flame  accelerations  will  be  looked  for. 

I 

PLANS  FOR  THE  COMING  YEAR 

’ The  systematic  study  of  initiation  using  the  CLOUD  program  will  con- 

tinue. After  a rather  complete  understanding  of  methane  initiation  is  ob- 
tained  we  plan  to  look  at  other  fuels  by  putting  in  their  kinetic  rate  laws. 

In  the  experimental  program  we  plan  to  study  in  some  detail  the  ef- 
fect of  various  sized  and  shaped  obstacles  on  the  acceleration  process  of  a 
free  flame  propagating  "through  a layered  combustion  mixture  which  is  open  at 
the  top.  This  is  a strictly  two-dimensional  study  so  it  should  be  easily 
amenable  to  analysis  by  hydrocode  techniques. 
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A SYNOPSIS  OF  RESEARCH  ON 
DDT  OF  GRANULATED  PROPELLANT 


Work  has  been  underway  at  the  University  of  Illinois  these  past  three 
years  to  analyze  the  unsteady  reactive  two-phase  flow  associated  with  DDT 
(Deflagration  to  Detonation  Transition)  of  confined  granulated  solid  pro- 
pellant and  explosives.  The  modeling  efforts  to  date  which  have  been  pub- 
lished are  listed  below: 

1.  Beckstead,  M.  W. , Peterson,  N.  L.,  Pilcher,  D.  T. , Hopkins, 

B.  D. , and  Krier,  H. , "Convective  Combustion  Modeling  Applied 
to  Def lagration-to-Detonation  Transition  of  1MX",  Combustion 
and  Flame  30,  231-241  (1977). 

2*  Krier,  H.  and  Kezerle,  J.  A.,  "Modeling  of  DDT  in  Granulated 
Solid  Propellant",  AF0SR-TR-78-07 , October  1977. 

3.  Krier,  H.  and  Gokhale,  S.  S. , "Modeling  of  Convective  Mode 
Combustion  through  Granulated  Propellant  to  Predict  Detona- 
tion Transition",  AIAA  Journal  16,  177-183  (February  1978). 

4.  Krier,  H. , Gokhale,  S.  S. , and  Hoffman,  S.  J.,  "Unsteady 
Two-Phase  Flow  Analysis  Applied  to  Def lagration-to-Detonation 
Transition",  AIAA  Paper  78-1013,  presented  at  the  14th  AIAA/SAE 
Joint  Propulsion  Conference,  Las  Vegas,  NV  (July  1978). 

The  formulation  of  the  theoretical  models  presented  in  the  above  list 
indicates  the  continuing  development  of  the  analysis  which  can  predict  tran- 
sition to  detonation.  From  the  countless  number  of  computer  simulations 
carried  out  we  have  determined  that  there  are  constitutive  relations  which 
are  "rate-determining"  functions.  In  addition  these  relations  are  generally 
not  well  defined  for  unsteady  flows  at  high  pressure  and  high  solids  loading  - 
the  regime  of  most  interest  to  the  DDT  problems. 

The  key  constitutive  relations,  in  some  representative  order  of  their 


* 

Appendix  A is  a technical  paper  extracted  from  this  report. 
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importance  (as  they  affect  the  predictions)  are: 

(1)  Gas-particle  viscous  interaction,  which  basically  determines  the 
hot-gas  permeability  into  the  unignited  portions  of  the  granu- 
lated bed. 

(2)  Gas-particle  heat  transfer  coefficient,  which  determines  the  con- 
vective energy  transfer  from  the  turbulent  flow  hot  gases  to 
particles. 

(3)  The  propellant  burning  rate  at  extreme  pressures  (of  the  order 
of  109  nt/m2)  and  rapid  rates  of  pressure  change,  dp/dt. 

(4)  The  ignition  criteria  that  fixes  the  time  that  the  particles  are 
ignited  by  the  convective  processes.  Heat  flux  rates  often 
exceed  10*  watts/m2. 

(5)  The  intergranular  stress  that,  under  highly  transient  conditions, 
limits  the  (extreme)  particle  compaction. 

(6)  The  gas  equation-of-state.  It  has  been  determined  that  a con- 
stant co-volume  correction  to  the  ideal  e.o.s.  is  not  valid  at 
these  high  pressures. 

(7)  The  temperature  dependency  on  the  specific  heat,  gas  viscosity, 
and  gas  conductivity. 

(8)  A DDT  criterion. 

In  the  study  reported  in  Ref.  4 above,  we  have  noted  that  the  gas- 
particle  drag  interaction,  as  generally  used  in  the  DDT  models,  can  lead  to 
sizeable  particle  motion  and  extreme  compaction,  unless  one  can  provide  the 
appropriate  intergranular  resistance  to  the  compaction.  In  short,  it  was 
not  possible  to  properly  predict  the  convective  mode  combustion  dynamics  in 
long  granulated  beds,  using  the  same  gas-particle  drag  interaction  law  that 
seemed  to  work  for  the  shorter  beds.  And  as  discussed  in  Ref.  4,  the  data 
base  that  has  been  used  to  correlate  the  pressure  drop  in  packed  beds  has 
been  carried  out  both  at  steady-state  conditions  and  at  Reynold's  numbers 
several  orders  of  magnitude  less  than  that  needed  for  the  DDT  flow  processes. 
In  addition,  these  correlations  were  based  on  tests  only  with  inert  particles. 


The  same  criticism  must  obviously  be  applied  to  the  heat  transfer  co- 
efficient, as  correlated  at  low  pressures  and  relatively  low  velocities,  re- 
sulting in  only  moderate  Reynold's  number  ranges.  From  the  sensitivity 
studies  reported  in  Refs.  2-4,  a general  conclusion  can  be  stated  that  both 
the  gas-particle  friction  coefficient  and  heat  transfer  coefficient,  when 
extrapolated  to  the  Reynold's  number  and  low  porosities  needed  in  the  DDT 
flows,  are  too  large,  by  at  least  one  order  of  magnitude. 

Propellant  and  explosive  linear  burning  rates  are  generally  expressed 
as  functions  of  the  ambient  pressure,  through  a power  law  steady-state  cor- 
relation. Measurement  of  burning  rates  are  generally  never  carried  out  at 
pressures  greater  than  10,000  - 20,000  psi.  Yet  the  pressures  predicted  to 
occur  in  the  confined  granulated  bed  during  the  flamespreading  can  exceed 
250,000  psi  (1.72  (103)  nt/m2) . And  the  pressure  can  change  with  time  at 
rates  of  the  order  of  1000  psi/ysec  to  10,000  psi/ysec!  Thus  it  is  question- 
able whether  steady-state  burning  rates,  extrapolated  from  much  lower  pres- 
sures,represent  the  dynamic  burning  rates  during  these  flamespreading  pro- 
cesses. Also,  the  gas  velocities  relative  to  the  particles  can  vary  from  10 
to  1000  meters/second,  bringing  unknown  factors,  such  as  erosive  burning 


augmentation  into  the  burning  rate  expression. 

One  also  requires  an  ignition  criterion  from  the  stimulus  of  the  con- 
vective heat  transfer  from  the  hot  gaseous  combustion  products  as  they  are 
forced  through  the  unignited  regions.  It  was  already  mentioned  that  heat 
fluxes  ranging  from  100  to  1000  BTU/in2  sec  are  calculated  to  occur  during 
the  DDT  process.  Little,  if  any,  data  exist  at  these  extreme  flux  rates  that 
tie  into  either  a critical  propellant  surface  temperature  or  a cril  cal  ig- 
nition energy.  With  these  flux  rates,  ignition  delay  times  are  often  less 
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EXPERIMENTS  TO  OBTAIN  INFORMATION 
ON  DYNAMIC  BED  PERMEABILITY 

The  pressure  loss  of  the  gases  flowing  through  packed  beds  of  various 
sized  solid  particles  has  been  investigated  by  numerous  researchers,  both 
theoretically  and  experimentally.  The  work  of  Sabri  Ergun  [5]  has  generaly 
been  accepted  as  the  most  accurate  and  reliable  in  this  field.  However, 
Ergun's  findings  are  only  valid  for  Reynold's  number.  Re,  (based  on  particle 
diameter  and  superficial  gas  velocity)  of  less  than  3000.  More  recently,  Kuo 
and  Nydegger  [6]  have  measured  flow  resistances  with  small  particle  packed 
beds  for  a Reynold's  number  range  between  500  - 15,000,  with  porosity  aver- 
aging 0.38. 

Figure  1 is  a plot  of  the  drag  coefficient,  as  correlated  by  Ergun, 

namely, 

f = {1.75  Re  + 150(1-(|))> 

pg  <P 

where  <f>  = porosity  (percent  volume  of  voids).  The  regime  where  Ergun  cor- 
related his  relation  is  shown  in  the  figure  along  with  the  region  of  interest 
for  the  DDT  conditions.  It  is  obvious  that  one  would  be  extrapolating  well 
beyond  a reasonable  range  when  using  Ergun' s (or  even  Kuo/Nydegger' s)  relation 
for  the  gas-particle  drag  coefficient. 

Even  more  recently,  Robbins  and  Horst  [7]  have  measured  steady-state 
flow  resistances  in  packed  beds  at  much  higher  pressures  than  used  in  Re- 
ferences 5 and  6 and  up  to  values  of  Reynold' s number  of  10**.  They  indeed 
found  that  the  friction  factor  at  the  higher  Reynold's  number  (Re  ^ 10**)  was 
about  50%  of  the  extrapolated  Ergun' s value  [5]  and  85%  of  the  extrapolated 
Kuo/Nydegger  value  [6] . 
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Our  experimental  test  facility  was  already  designed,  but  not  fabri- 
cated, about  the  time  when  Robbins/Gough's  work  became  available  for  us  to 
review.  In  the  experiment  to  be  described  in  the  following  pages,  we  hope 
to  verify  the  results  of  Ref.  7 and  also  carry  out  nonsteady  flow  tests. 

The  first  test  rig  is  described  below.  We  have  only  a limited  pressure 
range  to  work  with  at  this  time,  as  will  be  discussed. 

Experimental  Apparatus 

Using  ideas  from  those  investigators  having  previously  carried  out 

’ 

flow  resistance  measurements  an  experimental  apparatus  was  designed  and  con- 
1 structed  (See  Fig.  2).  The  preliminary  tests  are  being  conducted  using  com- 

pressed air  from  a battery  of  storage  tanks.  The  maximum  obtainable  pressure 
from  the  storage  tanks  is  450  psig.  Bourdon  gage  1 (Fig-  2)  is  a Solfrunt 
4.5"  test  gage  with  a range  of  0 - 1000  psig  and  with  a ±0.25%  of  span  ac- 
curacy. This  gage  is  used  to  monitor  the  pressure  in  the  storage  tanks.  A 
Grove  Pressure  Reducing  Regulator,  Model  82-829  (720  psig  max.  inlet),  allows 
us  to  regulate  the  compressed  air  to  the  desired  test  pressure.  The  com- 
pressed air  enters  a plenum  tank  which  contains  three  straightening  grids. 

This  is  designed  to  bring  the  gas  to  approximately  zero  velocity.  Bourdon 
gage  2 is  the  same  type  as  Bourdon  gage  1.  Gage  2 indicates  the  test  pres- 
sure, which  is  used  in  the  calculation  of  mass  flow.  At  the  end  of  the  plenum 
tank  is  a Sonic  Flow  Nozzle  from  American  Meter  Division,  Singer  Corporation. 
This  nozzle  has  a throat  diameter  of  0.375  inches  (and  an  accuracy  of  ±0.15% 
at  10  psig).  Opposite  the  Bourdon  gage  2 ic  an  Omega  Thermocouple.  It  is  a 
chromel  Alumel,  CASS-18G-12,  grounded  junction  probe,  that  provides  fast 
response  under  high  pressures.  The  Sonic  nozzle  opens  into  a straightening 
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section  that  contains  two  straightening  grids. 

The  test  section  has  a screen  mesh  at  the  entrance  and  a metal 
stainless  steel  grid  at  the  exit,  to  hold  the  packed  bed.  There  are  three 
rapid  response  pressure  transducers  mounted  on  the  test  section,  to  measure 
the  pressure  drop  of  the  gas  as  it  flows  through  the  packed  bed.  The  trans- 
ducers are  Setra  Systems',  Model  205,  with  a pressure  range  0 - 500  psig,  a 
one  millisecond  response  time,  and  an  accuracy  of  0.11%  at  full  scale. 

The  Omega  thermocouple,  as  described  previously,  monitors  the  tem- 
perature in  the  test  section  and  after  the  gas  has  exited  the  test  sections. 
Figure  3 gives  a more  detailed  description  of  the  experimental  apparatus 
from  the  plenum  tank  through  the  test  section. 

The  data  from  the  Setro  transducers  is  recorded  on  a Tektronix  5100 
series  Storage  Oscilloscope,  and  the  displays  are  permanently  recorded  with 
a Tektronix  C-5  Oscilloscope  Camera. 

The  readings  from  the  Omega  Thermocouples  are  individually  displayed 
on  an  Omega  Digital  Thermometer,  Model  2160A,  with  a response  time  of  less 
than  2 seconds,  and  a maximum  error  including  NBS  conformity  of  ±1°F.  These 
readings  and  the  Bourdon  gage  readings  are  recorded  in  a lab  book,  and  are 
used  in  conjunction  with  the  Sonic  flow  nozzle  to  calculate  the  mass  flow 
for  each  test. 


Limitations  of  Experimental  Measurements 


The  experimental  data  obtainable  from  this  apparatus  is  currently 
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limited  by  the  compressed  air  source  of  450  psig,  which  in  reality  gives  us  a 
workable  test  pressure  of  400  psig.  This  limitation  on  pressure  limits  the 
flow  rate  and  consequently  the  Reynold's  number  range  that  can  presently  be 
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studied.  It  is  not  possible  to  obtain  heated,  high  pressure  air  with  this 
apparatus.  Hence  we  are  not  now  able  to  study  the  heat  transfer  and  pres- 
sure drop  of  a hot,  high  pressure  gas  which,  in  fact,  exists  in  propellant 
beds. 

Ex-perimental  Procedures  and  Planned  Testing 

Tests  will  be  carried  out  initially  with  3 nun,  solid  glass  beads 
-3  3 

(density  2.73(10  ) g/mni  ) at  four  input  pressures,  namely  100,  200,  300 

and  400  psig.  Then  the  test  section  will  be  repacked  with  6 mm,  glass 

_3 

beads  (density  2.587  (10  ) and  flow  conditions  at  the  same  four  input  pres- 
sures. The  collected  data  will  be  analyzed  so  that  a valid  correlation 
between  the  coefficient  of  drag,  f , and  the  Reynold's  number,  and  hence 
the  mass  flow,  density  and  porosity  can  be  derived. 

Lead  shot,  the  same  size  as  the  3 mm  and  6 mm  beads,  but  of  different 
density,  will  be  also  used  in  the  test  section.  The  effects,  if  any,  of 
these  different  density  spheres  will  be  studied,  enabling  us  to  better  un- 
derstand the  parameters  that  affect  the  bed  permeability. 

During  these  tests,  the  pressure  transducers  will  initially  be  mounted 
vertically,  as  shown  in  Fig.  3.  Then  the  whole  test  section  will  be  rotated 
first  by  90°  and  then  180°  (so  that  they  will  be  located  on  the  bottom  of 
the  test  section) . This  will  enable  us  to  analyze  the  flow  under  identical 
test  conditions  from  three  axial  locations,  which  will  help  us  ensure  that 
the  flow  through  the  test  section  and  the  packing  of  the  beads  is  uniform. 

In  the  second  phase  of  this  experiment  we  will  pack  the  bed  with  inert 
cylindrical  propellant,  0.4026  inches  long  and  0.18332  inches  in  diameter 
and  repeat  tests  at  the  same  pressures  used  for  the  spherical  particles. 
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This  latter  study  will  attempt  to  verify  the  observations  made  by  Robbins 
and  Gough  [7]  on  the  geometrical  scaling  of  particles  on  the  bed  re- 
sistance. 

Our  planned  third  phase  is  to  mix  the  3 mm  beads  and  the  inert 
cylindrical  particles  and  then  record  the  pressure  drop  in  the  bed  versus 
the  percent  of  spherical/cylindrical  mixing,  to  determine  any  correlation 
between  the  pressure  drop  with  the  porosity  and  percent  size  mixing. 

The  fourth  proposed  phase  is  a transient  experiment.  This  will  con- 
sist of  attaching  several  gas,  high  pressure,  air  bottles  one  inch  down 
stream  of  the  regulator  and  allowing  the  gas  bottles  to  completely  empty 
from  2000  psig  to  about  15  psig.  The  recording  of  the  transient  pressure 
at  each  of  the  three  transducers  would  indicate  nonsteady  effects.  A plot 
of  this  transient  pressure  drop,  Ap,  should  cross  the  plot  of  the  steady- 
state  Ap  at  the  various  test  pressures  of  400,  300,  200  and  100  psig.  This 
will  enable  us  to  check  the  validity  of  our  coefficient  of  drag  correlation 
for  the  unsteady  flow  case,  which  is  what  actually  occurs  in  packed  beds  of 
interest  in  the  DDT  problems. 

Preliminary  Results 

We  are  just  beginning  to  take  data  now.  Our  initial  results  from  the 
tests  conducted  for  Ap  versus  mass  flow  of  air  lie  within  the  data  pre- 
sented by  Robbins  and  Gough  [7]  and  are  shown  in  Fig.  4.  Note  that  our  6 mm 
beads  which  are  a little  smaller  than  1/4  inch  (or  4/16)  used  in  Reference 
7,  provide  conditions  between  the  3/16"  and  5/16"  bead  size  data  of  Robbins 
and  Gough. 
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IMPROVEMENTS  IN  NUMERICAL  SCHEMES 
TO  PREDICT  DDT 

Summarized  below  is  a brief  description  of  our  efforts  and  progress  in 
improving  the  numerical  integration  scheme  which  is  utilized  to  solve  the 
DDT  model. 

(1)  Instead  of  using  the  operator  form  for  the  governing  equations  we  are 
anticipating  a shock  to  develop  as  transition  from  deflagration  to  detona- 
tion occurs.  The  simplest  form  of  the  Rankine-Hugoniot  relations  uses  the 
conservation  variables,  namely,  p,  pU,  pE  that  satisfy  the  jump  conditions. 
Hence  the  usage  of  conservation  variables  is  anticipated  to  produce  better 
results. 

(2)  We  have  also  revised  the  manner  in  which  we  specify  the  boundary  con- 
ditions at  the  closed  end.  The  correctness  of  the  boundary  conditions  for 
such  cases  is  insured  by  satisfying  zero  gradients  at  the  end  walls. 

(3)  For  the  predictions  of  deflagration  to  detonation  transition,  we  have 
prescribed  initial  conditions  so  that  the  problem  is  sure  to  start  as  a 
weak  deflagration.  Uniform  pressure,  temperature  and  loading  density  will 
lead  to  the  zero  velocity  profile  which  will  be  consistent  with  the  govern- 
ing equations.  The  problem  is  started  with  a perturbation  in  the  particle 
energy  by  assuming  a certain  fraction  of  bed  being  ignited  at  initial  time. 

(4)  An  attempt  was  made  to  understand  sensitivity  of  the  constitutive  re- 
lations for  the  particle-particle  interaction  force  and  its  bearing  in  v.he 
particle  energy  equation.  This  is  very  important  in  view  of  the  fact  that 
we  assume  that  the  particle  phase  represents  a continuum.  There  is  only  a 
limited  data  base  that  is  available  to  verify  predictions  of  particle. 
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heating  to  the  compression  of  granular  beds  under  high  pressure. 

(5)  The  numerical  method  which  we  have  used  so  far  is  the  same  one  we  have 
reported  previously,  i.e.,  an  explicit,  two-step,  predictor-corrector 
MacCormack  method.  In  the  current  study  the  shock  is  not  dealt  with  any 
special  treatment.  This  method  can  be  termed  as  a shock  smearing  technique. 
Since  the  governing  equations  are  inviscid  in  nature,  it  is  necessary  to 
incorporate  damping  to  insure  numerical  stability.  Due  to  ‘•he  higher  non- 
linearity and  interdependent  phase  quantities,  we  found  it  necessary  to  use 
artificial  viscosity  similar  to  Von  Neumann  and  Richtmyer  type  all  the  time. 
After  trying  out  various  types  of  damping  functions  we  found  that  the  three 
point  weighted  average  is  relatively  easy  to  use  and  provides  quite  satis- 
factory predictions.  Figure  5 shows  the  pressure  distribution  at  various 
times  for  a typical  case.  It  shows  a gradual  steepening  of  the  pressure 
wave  which  could  lead  to  a strengthening  of  the  shock.  Figure  6 shows  loci 
of  various  fronts.  As  could  easily  be  seen,  the  flame  front  is  accelerating 
moderately  throughout  the  last  half  portion  of  the  propellant  bed. 

(6)  We  are  in  the  process  of  developing  another  differencing  analog.  Lax- 
Wendorff.  It  is  not  as  yet  adequately  tested  and  hence  the  results  are  not 
reported  in  this  memo. 

(7)  We  also  attempted  nondimensionalizing  all  of  the  governing  equations. 
Due  to  the  interdependency  of  the  two  phases  this  poses  great  difficulties 
and  only  limited  success  in  solution  of  the  equations. 

(8)  We  are  also  in  the  process  of  developing  difference  analog  for  multi- 
size particles.  In  the  finished  stage  our  DDT  computer  program  will  have 
the  capability  to  handle  a wider  variety  of  problems  with  the  propellant 
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particles  of  various  sizes,  as  well  as  greater  solids  loadings.  Even 
though  the  logic  of  developing  such  a program  is  fairly  straightforward, 
with  the  usage  of  two-dimensional  storage  arrays,  such  a code  poses  dif- 
ficulties. Most  notable  among  them  are  the  definitions  and  the  distri- 
bution of  the  interphase  transfer  quantities  such  as  drag  and  heat  transfer. 
There  is  also  a problem  in  defining  the  ignition  front  through  multisize 
reactive  particles.  Additional  insights  are  needed  in  order  to  make  this 
analog  successful.  Figures  7 and  8 demonstrate  the  results  of  our  initial 
efforts  in  developing  the  multisize  particle  program.  Figure  7 indicates 
that  much  higher  pressures  occur  at  the  corresponding  times  for  the  multi- 
size particles  as  compared  to  unisize  particles.  Figure  8 shows  the  locus 
of  flame  front  for  the  two  cases. 

(9)  Ir.  the  past  our  attempts  were  restricted  to  fixed  space  Eulerian  mesh. 
During  this  research  year  we  intend  to  develop  a Lagrangian  formulation 
which  would  allow  us  rezoning  and  restructuring  the  mesh  depending  on  the 
strength  of  the  shock.  If  this  technique  is  successful,  we  should  be  able 
to  clearly  predict  a detonation  which  follows  from  an  accelerating  defla- 
gration through  a porous  reactive  bed  of  propellant. 
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SUMMARY 

A two-phase  reactive  flow  model  is  derived  and  solutions  are  applied 
to  the  analysis  of  the  unsteady  convective  heat  transfer  in  initially  packed 
beds  of  granulated  solid  propellant  or  explosives.  The  resulting  pressure 
wave  and  flame  spreading  event  is  considered  the  physical  process  that  may 
provide  for  deflagration-to-detonation  transition  (DDT)  in  such  granulated 
beds.  Discussions  are  included  as  to  the  appropriate  constitutive  laws  re- 
quired to  complete  the  hydrodynamic  model.  The  paper  concludes  with  a brief 
assessment  regarding  both  the  limitations  of  the  model  and  another  possible 
DDT  mechanism. 
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INTRODUCTION 

The  transition  from  combustion  to  detonation  in  condensed,  but  porous 
propellants  and  explosives  has  in  recent  years  received  attention  both  from 
analytical  and  experimental  efforts.  A large  amount  of  the  experimental 
work  has  been  carried  out  by  Soviet  scientists,  much  of  which  has  been  docu- 
mented by  A.  F.  Belyaev  et_  al.  [1]  and  recently  reviewed  by  H.  H.  Bradley 
and  T.  1.  Boggs  [2J.  The  work  of  R.  R.  Bemecker  and  D.  Price  13,4]  repre- 
sents detailed  and  specific  experimental  research  toward  understanding  de- 
flagration-to-detor.ation  transition  (DDT)  in  porous  explosives.  These 
latter  studies  clearly  indicate  that  the  buildup  to  detonation  represents  a 
coupling  between  the  pressure  (shock  or  compression)  fronts  and  the  con- 
yectively  driven  flame  front. 

As  would  be  expected  if  detonation  is  to  follow  ignition  (and  defla- 
gration) it  is  necessary,  above  all,  to  ensure  a sufficiently  rapid  pressure 
buildup.  In  a porous  combustible  system,  this  is  achieved  as  a result  of 
the  penetration  of  the  gaseous  combustion  products  into  the  porous  interior 
which  leads  to  accelerating  ignition  of  additional  burning  surfaces.  Thus 
heat  transfer  by  conduction  is  replaced  by  convective  heat  transfer.  The 
specific  mechanism  to  achieve  DDT  is  yet  to  be  fully  understood,  but  pro- 
posed phenomenological  mechanisms  have  been  proposed  by  Bemecker  and  Price 
[3] , by  R.  W.  Van  Bolah  et  al  |5] , and  by  Belyaev  et^  al  [1] . 

The  analytical  work  to  describe  the  DDT  process  by  solving  the  full 
reactive  hydrodynamic  conservation  equations  has  been  underway  at  the  Uni- 
veristy  of  Illinois  for  the  past  few  years  and  is  reported  in  Refs.  6-9.  The 
model  by  H.  Krier  and  S.  S.  Gokhale  [8]  centered  upon  a two-phase  continuum 
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knowledge  of  the  reaction  kinetics  in  the  solid, formation  of  active  (hot) 
sites,  and  the  appropriate  equation  of  state  for  expressing  the  solid  density 
in  terms  of  the  extremely  high  pressures  and  temperatures.  This  work  does 
not  include  these  stages  of  the  transition  process,  but  our  calculations 
show  that  with  a rapid  transition  of  the  convectively  driven  flame  front, 
in  conjunction  with  the  interaction  of  the  pressure  fronts,  one  may  approxi- 
mately indicate  conditions  which  might  lead  to  detonation. 

THE  MODEL:  FIELD  BALANCE  CONSERVATION  EQUATIONS 

The  one-dimensional  conservation  equations  for  a two-phase  reactive 
mixture  have  been  derived  by  most  investigators  through  the  concept  of 
separated  flow  as  defined  in  the  text  by  G.  B.  Wallis  [15].  This  approach 
considers  two  distinct  flows,  each  through  a separate  control  volume,  such 
that  the  sum  of  the  volumes  represents  an  average  mixture  volume  and  the 
sums  of  the  properties  of  each  flow  represent  average  mixture  properties  of 
the  fluid.  With  this  approach,  separate  equations  for  continuity,  momentum, 
and  energy  are  written 'for  each  phase  and  are  solved  with  a state  equation  to 
describe  the  overall  flow.  Panton  [16]  has  provided  the  rigorous  formula- 
tion with  the  important  assumptions  required  to  express  these  field  balance 
equations.  There  are,  of  course,  many  other  analyses  that  have  arrived  at 
similar  forms  of  the  equations.  The  Appendix  presents  an  outline  of  the 
method  we  have  used  to  work  out  our  conservation  equations. 

Assumptions 

The  important  assumptions  upon  which  this  separated  flow  model  is 
based  are  (1)  The  two  phases  are  assumed  interdispersed , but  separate, 
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coupled  only  by  the  appropriate  interaction  terms;  (2)  Each  phase  is  a con- 
tinuum and  therefore  derivatives  are  uniquely  defined;  (3)  The  total  cross 
sectional  area  is  the  sum  of  the  gas  and  solid  flow  areas.  (The  problem  is 
quasi-one-dimensional . ) ; (4)  Solid  particles  are  large  compared  to  molecules 
and  hence  do  not  contribute  to  the  total  pressure  of  the  mixture;  (5)  When 
combustion  of  particles  causes  mass  transfer 'between  the  phases,  the  solid 
phase  always  loses  mass  while  the  gas  phase  gains  it,  i.e. , T ^ 0;  (6)  The 

diameter  of  a typical  particle  can  be  thought  of  as  the  average  diameter  of 
the  total  number  of  particles;  (7)  All  gases  obey  the  nonideal  Noble-Abel 
equation  of  state  with  a variable  co-volume;  (8)  Heat  transfer  to  the  par- 
ticles due  to  conduction  and  radiation  is  neglected;  (9)  The  density  of  the 
* 

solid  phase  is  constant;  (10)  The  equations  are  laminar  in  the  sense  that 

t 

turbulent  flow  due  to  the  two-phase  nature  of  the  fluid  has  been  averaged 
out.  (See  Panton  {16]  for  additional  discussions  regarding  this  assumption.); 
(11)  The  gases  are  inviscid  except  for  their  action  on  the  particles  through 
the  drag  tern;  (12)  The  fluid  properties  of  and  Cv  are  assumed  constant, 
but  gas  viscosity  and  conductivity  are  generally  known  functions  of  tem- 
perature. „ 

The  Conservation  Field  Balance  Equations 

Some  of  the  details  in  the  drivation  of  the  six  field  balance  equa- 
tions have  been  presented  in  the  Appendix.  In  order  to  specifically  describe 
the  state  of  a two-phase  fluid  one  needs  equations  to  solve  for  the  following 

variables:  p,u,T,$,u,  and  T . The  state  equations  are  added  to  the 
Egg  P P 

field  balance  equations  to  obtain  the  seven  required  equations.  Therefore 
two  continuity  equations,  two  momentum  equations,  and  two  energy  equations 
will  be  necessary.  These  equations  are  listed  below  in  operator  form. 
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It  has  been  duly  observed  that  several  investigators,  in  particular 
Kuo  et  al.  [10,11]  and  Gough  [12],  do  not  write  a solid  phase  energy  balance 
symmetric  with  the  gas  phase  energy  balance,  like  the  Eq.  (6)  above.  In- 
stead, these  investigators  have  decided  to  simply  solve  the  unsteady  heat 
conduction  equation  for  the  solid  particle  and  monitor  the  heat  gained  from 


6 


L 


the  gas  in  the  fom  of  convective  heat  transfer.  In  this  way  particle  igni- 
tion can  be  tied  to  a critical  surface-ignition  temperature..  That  is,  one 
solves  for  the  temperature  distribution  within  the  particles. 

Since,  in  the  spirit  of  the  one-dimensional  analysis,  we  wish  to  re- 


tain only  one  variable  to  define  particle  internal  energy,  i.e.,  C T , our 

VP  P 


approach  requires  retention  of  a particle  energy  equation,  e.g.,  Eq.  (6a). 
One  could  alternately  derive  the  particle  energy  equation  by  a simple  energy 
balance  on  one  particle,  and  then  correct  the  relation  to  account  for  a 
(!-<*>)  solids  fraction.  In  that  case  one  gets. 


Solid  Energy 
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Note  that  Eq.  (6b)  may  actually  be  more  appropriate  than  Eq.  (6a),  since 
the  form  of  Eq.  (6a)  implies  that  an  acceleration  of  the  particles  decreases 
the  particle  energy,  that  is,  the  solid  phase  is  a pseudo-fluid. 

In  Eqs . (1-6) , 
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The  subscript  "g"  denotes  gas  and  "p"  denotes  particle.  The 


mixture  cross-sectional  area  is  denoted  by  S,  where  S = A + A . 

P g 


Also  the  gas  generation  term  for  spherical  particles  of  radius  rp  is  noted  as 


T = - T where  , 
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and  the  propellant  burning  rate,  r,  is  assumed  given  by  the  relation 


T 

m(  1 

_P_ 

_E_| 

T 

p 

lPoj 

°J 

where  t^,  a,  and  n are  assumed  known  constant. 

Actually  propellant  burning  rates  at  the  high  pressures  required  here  have 
generally  never  been  measured.  Dynamic  burning  rates,  i.e. , r = r(dp/dt)  may 
also  be  important,  but  there  are  no  data  to  take  this  into  account. 

The  equation  of  state  for  the  propellant  gases  at  the  expected  high 


pressures  is  given  by 


P[-  - B(p) ] = RT 


where  B(p)  is  variable  co-volume  correction  based  on  the  work  of 
M.  Cook  [17]. 

The  interphase  heat-transfer  $ represents 

$ = J.  (i_(t))h  (T  - T ) 
rp  Pg  g P 

and  h = heat  transfer  coefficient  based  on  Denton's  work  [18]  through 
packed  beds  of  inert  particles.  We  therefore  assume 


h = 0.58  Re0,7Pr'33  (10) 

Pg  rp 

Likewise  a constitutive  relation  is  required  for  the  interphase 
viscous  forces,  given  by  the  term  F in  Eqs.  (3)  and  (4).  Here  we  have 
utilized  either  Ergun's  drag  correlation  [19]  or  Kuo/Nydegger 's  relation  [20], 
That  is 

7 ’ **  (Ug  - V <lla) 
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where  for  Ergun's  data. 


and  the  Reynold  number, 
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It  should  be  noted  that  both  the  coefficients  h and  f were  correlated 

Pg  Pg 

at  pressures  and  temperatures,  and  thus  at  Reynolds  numbers,  several  orders 
of  magnitude  less  than  for  conditions  arrived  at  in  the  DDT  problem.  In 
addition,  the  solid  particles  here,  once  ignited,  generate  gases  at  their 
surface,  so  that  both  Denton's  and  Ergun's  correlations  may  only  marginally 
apply. 

Additional  relations  are  required  for  the  gas  conductivity,  k , gas 

g 

viscosity,  y , and  gas  specific  heat,  C , as  they  vary  with  temperature 

8 Vg 

and  pressure,  and  are  assumed  known. 

Finally,  a relation  is  required  for  the  resistance  to  compaction,  at 

high  solid  loading,  as  symbolized  by  the  term  T^.  Again  very  little  data 

are  known  for  the  dynamic  normal  axial  stress  for  an  aggregate  of  particles, 

given  by 

but  we  have  used  the  same  expression  / Kuo  and  Summerfield  [10]  have  utilized. 

That  is,  beyond  some  critical  solids  load,  less  than  (l-<j>  ), 

c 


(12) 


1-4.  / 1-4>C  l-4> 


where  K is  a bulk  modulus  for  the  aggregate.  Since  for  uni-sized 
spherical  particles,  minimum  compaction  cannot  exceed  74%,  the  value  for  K 
must  be  large  enough  so  that  the  normal  axial  stress,  proportional  to  St^/Sx, 
prevents  further  acceleration  of  the  particles. 
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Problem  of  Concern 


The  two-phase  separated  flow  model  generates  a set  of. six  nonlinear, 
inhomogeneous,  and  coupled  hyperbolic  partial  differential  equations,  as 
demonstrated  above.  The  numerical  method  chosen  for  solving  these  equations 
is  the  explicit  two-step  MacCormack  scheme.  This  method  was  chosen  for  its 
accuracy  and  because  it  is  well  documented  as  a two-step,  predictor-corrector 
method  of  integration.  Ref.  9 gives  details  of  how  we  utilized  this  scheme, 
how  stability  was  assured,  and  how  the  boundary  conditions  were  handled  at 
the  closed  ends.  Equations  (1-6)  have  also  been  shown  to  constitute  a 
hyperbolic  system  of  equations,  as  required. 

Typically  a bed  of  highly  loaded,  energetic  small  particle  propellant 
with  solids  loading  ranging  from  0.5  to  0.7  was  considered  and  analyzed  with 
closed  ends.  Bed  lengths  from  5 cm  to  15  cm  were  considered,  and  the  space 
domain  divided  in  grids  of  ix  = 0.1  cm.  Stability,  dependent  upon  the  gas 
velocities  and  mixture  sound  speed  required  that  the  time  increment.  At,  be 
of  the  order  of  0.5  y seconds. 

The  problem  was  initialized  by  assuming  that  at  time,  t = 0,  a pres- 
sure pulse,  typically  10  M N/m2  (1500  psi)  was  distributed  over  the  first 
10  percent  of  the  bed.  It  was  also  assumed  that  near  one  end  3%  of  the  bed 
was  initially  ignited,  that  is  T^  > Tp^gn)'  The  ignition  criteria  employed 
here  requires  that  if  the  particle  energy  > ®p(ign)*  an  assume<*  hnown 
value,  they  begin  to  burn,  at  a rate  given  by  r^.  Thus,  - Ep(ign)/Cy 

After  several  hundred  independent  calculations  a baseline  case  with 
standard  input  conditions  was  deduced.  It  should  be  noted  that  the  solid 
phase  stress  tensor  in  Eq.  (4)  does  not  include  the  gas  pressure.  The  reason 
for  this  is  that  it  has  been  shown  [8]  that  the  set  of  equations  will  be 
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totally  hyperbolic  only  if  there  is  no  gas  pressure  gradient  in  the  solid 
phase  momentum  equation. 

Table  1 summarizes  the  most  important  input  used  to  carry  out  the 
typical  results  shown  below.  The  main  results  presented  from  this  analysis 
are  the  pressure  distributions  as  time  progresses,  the  flame  front  locus 
and  pressure  front (s),  as  well  as  developments , in  the  temperature  and 
velocity  fields. 

Results  from  the  Baseline  Case 

Figures  1 through  5 present  the  distribution  histories  of  the  fluid 
dynamic  variables  in  a bed  of  solid  propellant  7.56  cm  (3  in.)  long,  with  a 
solids  loading  at  60  percent  (<j>o  = 0.40).  The  results  in  these  figures 
were  termed  the  baseline  case,  as  defined  above. 

Figure  1 outlines  the  calculated  pressure  distributions  at  four  dif- 
ferent times.  As  in  previous  work  related  to  this  problem  (see  Ref.  8) 
the  appearance  of  a "continental-divide"  in  the  interior  of  the  bed  is  pro- 
nounced and  pressure  gradients  are  extremely  steep.  In  addition,  the  average 
pressures  behind  the  front  are  very  high,  reaching  in  excess  of  3.4  G N/m2 
(0.5  x 106  psia).  Such  pressures  may  be  unrealistically  high  when  compared 
to  experiments  in  closed  pipes  as  reported  in  Refs.  3 and  4,  due  to  the  fact 
that  here  we  are  assuming  incompressible  solid  particles,  all  spherical  in 
shape,  so  that  once  ignited  all  surfaces  generate  gases.  Irregular  shaped, 
compressible  particles  may  compact  in  such  a way  to  prevent  all  available 
surfaces  from  burning,  even  though  the  solid  temperatures  are  beyond  the 
critical  ignition  temperature. 

The  locus  of  the  flame  (ignition)  front  for  the  baseline  case  is  pre- 
sented in  Fig.  2,  where  the  deflagration  speed  (slope  of  the  x-t  locus) 


accelerates  to  several  nun/ps.  Also  shown  in  this  figure  is  a line  labeled 
the  "pressure  front."  This  locus  was  derived  by  noting  the  midpoint  of 
the  pressure  wave  shown  in  Fig.  1 at  various  times.  It  appears  that  the 
pressure  front  and  the  flame  front  are  coincident  as  the 


11 


Table  1 

TYPICAL  INPUT  DATA 


Parameter 


Value 


Initial  bed  temperature,  T^,  Tp 
Bulk  itnition  temperature,  T^n 
Range  of  initial  bed  porosity,  <t> 

Propellant  burning  rate  proportionality  constant,  b 
Propellant  burning  rate  index,  n and  m 
Propellant  density, 

Initial  grain  diameter,  d = (2r  ) 

R R C Vi  pm 

Chemical  energy  released,  E , = (E  -E  ) 

chera  g p 

Latent  heat  of  propellant  sublimation  (E  ) 

P 

Molecular  xceight  of  gas,  MW 

Covolume  of  propellant  gas,  (at  moderate  pressures) 
Specific  heat  ratio  of  gas,  y 
Gas  viscosity,  yg  (at  2000  K) 

Gas  specific  heat  at  constant  volume  (2000  K)C 

a 

Solid  specific  heat  at  constant  volume,  C 
Total  bed  length,  ** 

Bulk  modulus,  K (Eq.  T.2) 

Compaction  initiation  porosity,  <J>c  (Eq.  12) 


294. 0°K 
314  °K 

0.30  - 0.50 
1.24  cm/s  (MPa)n 
0.90  and  0.25 
1.58  g/cm3 
400  ym 
5.48  MJ/kg 
0.22  MJ/kg 
22.6  kg/kg-mole 
1.08  cm3/g 
1.252 

4.45  N-s/m2 
0.340  cal/gram  °K 
0.305  cal/gram  °K 
75.0  mm 
48  M N/m2 
0.42 
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flame  accelerates  to  complete  Ignition. 

Figure  3 shows  both  the  gas  velocity  and  the  propellant  particle 
velocity  distributions  at  the  chosen  times.  As  might  be  expected,  the 
velocity  peaks  shown  are  a consequence  of  the  steep  pressure  front  arriving 
at  the  various  locations.  Also  worth  noting  are  the  extreme  gas  velocities, 
often  exceeding  1500  a/s,  and  the  fact  that  peaks  in  particle  velocity 
lag  peaks  in  gas  velocity  at  any  given  time. 

Figure  4 presents  gas  temperature  and  particle  temperature  distribu- 
tion histories.  With  the  chemical  energy  fixed  for  this  propellant,  the  gas 
temperatures  are  generated  at  values  of  roughly  3600°K.  These  temperatures 
decay  as  the  gases  lose  heat  to  the  solid  by  convection  into  the  bed  interior. 
Later,  as  the  pressure  front  steepens,  the  gases  are  compressed  in  the  in- 
terior to  values  far  exceeding  their  incoming  values.  The  subsequent  par- 
ticle temperatures  are  increased  not  only  due  to  the  convective  heat  transfer, 
but  also  due  to  the  rapid  changes  in  kinetic  energy  of  these  mobile  particles? 
The  last  predictions  displayed  for  the  baseline  case  are  the  porosity  dis- 
tribution histories  shown  in  Fig.  5.  The  porosity  profiles  show  that  neav 
the  front  of  the  bed„(which  is  partially  ignited  at  t ■ 0)  the  porosity  in- 
creases fairly  rapidly.  This  is  due  to  both  a reduction  in  volume  of  the 
particles  due  to  burning  and  the  forward  motion  of  these  particles  being 
dragged  along  by  the  gases.  Subsequent  to  this,  the  particles  are  compacted 
somewhere  in  the  bed  interior  to  solids  loadings  greater  than  the  original 
60  percent.  (Recall  that  the  solid  fraction  equals  (1~4>).)  Also,  for  the 
results  shown  here  a significant  fraction  of  the  bed  has  been  consumed  (due 
to  rapid  burning  at  the  high  pressures)  at  times  in  excess  of  50  Us. 

Decreasing  the  particle  size  from  the  200  pm  diameter  (r^  = .004  inches) 

With  Eq.  (6b)  used  instead  of  Eq.  (6a),  particle  energies  are  latered  only 
by  convective  heat  transfer  to  the  particle  surface. 


yields  a trend  similar  to  increasing  the  solids  loading.  In  the  case  shown 

in  Fig.  6,  decreasing  the  radius  to  100  pm  results  in  an  increase  in  the  (S/V) 

ratio  of  the  particles  of  100  percent,  since  (S  /V  ) = 3(l-«j>  )/r  . For  the 

p Cl  o op 

input  data  base  used  here,  the  predicted  pressures  for  the  small  particles 
size  exceeded  7*(109)  N/m2.  Therefore  the  maximum  pressure  in  the  bed  as  a 
function  of  particle  size,  as  presented  in  the  insert,  is  shown  for  an  ar- 
bitrary time  of  26.6  ys,  and  not  when  the  whole  bed  was  ignited. 

An  alternate  pressure  wave  "front"  can  be  defined  by  first  plotting 
pressure  vs.  time  at  any  of  the  downstream  x-locations  in  the  bed,  and  this 
is  shown  in  Fig.  7 for  a case  using  a slightly  different  form  of  the  par- 
ticle energy  equation,  i.e. , Eq.  (6b).  The  locus  of  a pressure  front  can 
then  alternately  be  defined  as  the  time  at  any  position  vrhen  the  pressure 
begins  its  rapid  increase.  For  example,  on  Fig.  7,  at  x = 3.8  cm,  the  pres- 
sure front  is  said  to  arrive  at  82  Jis.  This  locus  is  plotted  along  with  the 
ignition  (flame)  front  on  Fig.  8.  Notice  the  fairly  abrupt  change  in  flame 
(deflagration)  spead  at  approximately  60  ys,  which  approximately  coincides 
with  the  intersection  of  the  pressure  front  as  obtained  from  Fig.  7:  The 
calculations  of  the  £larae  front,  pressure  wave  front  shown  in  Fig.  8 have 
many  of  the  same  features  as  those  observed  in  the  DDT  tests  reported  by 
Bernecker  and  Price  13] . Although  the  predictions  shown  here  cannot  yet  be 
termed  actual  DDT,  they  do  represent  indications  of  the  rapid  propellant 
combustion,  producing  high  pressures  which  drive  hot  product  gases  into  the 
unignited  pores  in  a rapidly  accelerating  fashion.  When  calculations  of  the 
type  shown  in  the  previous  figures  were  repeated  for  lov-to-medium  energy 
content  solid  propellant,  say  750  cal/g,  no  such  dynamic  acceleration  (as 
shown  in  Fig.  8)  were  predicted,  even  in  bed  lengths  twice  the  nominal  value 
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of  7.5  cm.  Similarly  flame  speeds  with  the  high-energy  propellant  (see 
Table  1),  but  with  a burning  rate  constant,  b,  one-half  of  the  nominal 
baseline  value,  never  accelerated  beyond  0.2  mm/ys  in  the  bed  lengths 
studied  here. 

CONCLUDING  REMARKS 

Using  the  concept  of  a separated— flow  two-phase  mixture  of  solid 
propellant  particles  and  their  resulting  gaseous  products,  the  one-dimensional 
form  of  the  appropriate  governing  conservation  relations  were  solved  along 
with  a series  of  constitutive  relations.  It  was  mentioned  that  many  of  these 
latter  relations  were  utilized  for  conditions  in  this  problem  which  are  well 
removed  from  those  as  they  were  originally  obtained.  The  report  by  Krier 
and  Kezerle  J9]  presents  various  sensitivity  studies  in  which  the  inter- 
phase heat  transfer  and  drag  correlations  were  altered,  and  where  the  ignition 
energy  was  varied.  The  conclusions  derived  from  such  a study  were  that,  on 
the  whole,  the  results  shown  here  represent  the  general  behavior  of  the  flow 
variables  during  the  dynamic  flame-spreading  process.  This  does  not  mean 
that  one  is  necessarily  satisfied  with  the  accuracy  of  these  constitutive 
laws.  But  rather,  if  precise  predictions  for  DDT  are  to  be  reported,  ad- 
ditional two-phase  experiments  (at  high  pressures  and  high  gas  temperatures) 
will  first  be  required. 

The  key  conclusions  from  this  study  are  then, 

(1)  On  the  average,  the  convective  mode  combustion  dynamics  through  granu- 
lated propellant  is  the  same  using  either  the  separated  flow  concept  or 
the  continuum  mixture  concept  of  Ref.  8,  although  in  this  paper  we 


have  not  presented  such  comparisons. 

(2)  As  would  be  expected,  only  highly  energetic,  rapid  burning  rate 
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propellants  exhibit  the  acceleration  to  detonation-like  flame  and 
pressure  fronts. 

(3)  Quantitative  predictions  for  the  run-up  distances  to  detonation  may 
have  to  be  delayed  until  better  information  is  available  for  the 
gas-particle  forces  and  heat  transfer,  for  the  particle-particle 
forces,  and  for  the  solid  ignition  criteria  for  these  transient  flow 
processes. 

(4)  However,  with  pressures  in  the  granulated  bed,  typically  predicted  to 
exceen  20  Kbar  (.2  G N/m2)  , in  time  periods  under  50  ys,  one  may  begin 
to  consider  DDT  phenomena  of  the  condensed  phase  alone,  due  to  clas- 
sical shock  initiation.  This  aspect  was  not  considered  in  this  study. 
If  this  were  a possibility,  one  might  consider  the  existence  of  an 
alternative  hazard  potential  for  solid  propellants  or  explosives,  in 
which  a portion  of  the  solid  first  sustains  damage  in  which  cracks 
formed  which  were  partially  filled  with  small  granulated  particles. 

If  such  a particle  laden  fissure  finds  a source  of  ignition,  then  the 
convective  mode  combustion  of  these  small  particles  under  confinement 
can  generate  an^accelerating  flame  and  the  subsequent  high  "shock" 
pressures  to  detonate  the  compacted  consolidated  solid  particles  in  the 
crack.  In  order  to  analyze  such  a shock  hydrodynamic  condition  one 
would  couple  the  solution  of  the  granular  bed  flame  spreading  with  the 
dynamic  deformation  of  the  compacted  solid.  An  equation  of  state  would 


be  required  for  the  solid  density. 


p , as  a function  of 
P 


the  pressure 


and  temperature.  Then  a shock  analysis  of  transition  to  detonation 


would  follow,  as,  for  example,  by  the  formalism  outlined  by  Forest  [22]. 
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APPENDIX 


Some  Details  of  the  Derivation  of  the  Field-Balance 
Conservation  Equations  for  Two-Phase  Separated  Flow 


There  are  many  ways  to  derive  the  quasi-one-dimensional  continuity. 


momentum,  and  energy  equations  for  the  gas  and  particle  phase.  The  work  of 
Panton  J[16J  is  a good  example  of  a rigorously  mathematical  approach.  An 


alternative  but  closely  related  derivation  is  shown  below  in  which  the  gas 
phase  and  solid  phase  are  treated  separately,  in  their  own  control  volumes, 
such  that  the  sum  of  these  volumes  adds  to  the  mixture  volume.  But  each 


phase  is  treated  as  a continuum,  so  that  when  the  solid-phase  is  analyzed 


it  is  treated  as  a "fluid"  of  properties  p , u , and  T 


P P 


Conservation  of  Mass: 


For  the  gas  phase  one  begins  by  expressing  gas-continuity  as 


o(p  A ) 9(p  u A ) 


+ r S 

c 


(A-l) 


where  A = 6S  and  S is  the  total  cross-sectional  area  of  the 
g 

(mixture)  control  volume.  Defining  p^  = 4>P  » one 


3pi  1 3ua  9p 

IT  + S Pl“g  + P1S  3^  + U8S 


(A-2) 


The  particle  phase  is  treated  in  a similar  manner,  except  u replaces  u , 

P © 

A is  replaced  by  A , where  A = (1-40 S and  T = - T . Here  T is  the  gas 
g P P CP  C C 

source  term  due  to  pyrolysis  of  the  solid  particles.  One  arrives  at 


IS 


-l  r\q 

“ + S p2up  3^  + p2s  3^  + ups  dT 


(A- 3) 


Conservation  of  Momentum 


where  p„  = (l-<j))p  . 

2 P 


For  the  gas  phase  one  simply  expresses  the  time  rate  of  increase 
of  momentum  inside  the  control  volume  to  be  equal  to  the  sum  of  all  forces 
acting  on  the  control  surfaces.  Thus, 


3(p  A u ) 3(p  A u2) 


Op  

+ ruS-A  -f^-FS 
c p g 3x 


(A-4) 


Note  that  it  is  assumed  the  pressure  gradient  force  acts  only  on  the  area 
A . It  is  also  assumed  that  since  the  gas  is  formed  at  a velocity  equal  to 

O 

the  particle  velocity,  there  is  a source  of  gas  momentum  per  unit  area  equal 
to  ^cUp’  0*  course  the  solid  phase  must  have  a sink  of  momentum  equal  to 
that  amount.  Also,  here,  F represents  the  viscous  interaction  forces  be- 
tween the  gas  and  all  the  solid  particles.  j 

To  express  (A-4)  in  operator  form,  one  again  introduces  the  variable 

p,  then  multiplies  Eq.  (A-2)  by  u , then  subtracts  that  equation  from  (A-4) 

S 


to  get 


3u  3u  » _ 

P, + P.u  = T (u  -u  ) - p v F 

1 ot  1 g 3x  c p g 3x 


(A-5) 


For  the  solid  phase  one  begins  with  the  momentum  balance  similar  in 
form  to  Eq.  (A-4),  i.e.. 


3IpC1-P)Su1  3[p  (1-<}>)Su2J 


- T u S - (l-p)S  -r-x-  + FS 
c p 3x 


(A-6) 


Again  using  the  Eq.  (A-3) , first  multiplied  by  Up  and  substituting  into 
(A-6)  one  arrives  at 
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3u  3u  3t  _ 

p2  d^+  P2Up  3^=  3^+F 


(A- 7) 


where  we  assume  p 14  constant,  and  where  T represents  the 
P P 

normal-axial  compressive  stress  (see  Eq.  12). 


Conservation  of  Energy 

For  the  gas  phase,  the  first  law  of  thermodynamics  for  a control 
yolume  proyides  that 


p A q = Fu  S + OS  + [p  A (E  +u2/2)]  + ~ 
Hg  g e p 3t  lhg  gv  g g 3x 


•)  > 
u2  P 

P A 

E +-£•  + — 6 

u 

g g 

g 2 P 

g 

{ 8j 

- r S 

c 


u5l 

chem  _£ 
?g  2_1 


(A-8) 


where 


E = fc  dl 
g v 


qe  “ 3x 
,chem 


l 3*J 


g 


(heat  conduction)  (neglected  later) 
= heat  of  combustion 


T S -j-  = kinetic  energy  added  because  gases  are  generated  with 
' the  velocity  u^. 


P = interphase  energy  transfer  from  gas  to  particles  and 
interphase  dissipation  work 


q _ + (T._ U_<f0 


pg  3x  xx  g 


, where  the  gas-shear  stress  T 


is  neglected  later 


xx 


and  where 


q = h (T  -T  )<J>(S  /V  ) 
Mpg  Pg  g P pm 


where  (S  /V  ) ■ 3 (1— <f>)  / r . 

pm  p 


To  arrange  Eq.  (A-8)  in  operator  form,  one  subtracts  from  it  the  product  of 


(E  + u2/2)  times  the  continuity  Eq.  (A-2),  and  the  product  of  (u  ) times 
g g ® 


r 1 l IT  i m*'  - -• 
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the  momentum  Eq.  (A-5)  to  get 


u n~— ' 
g 3x 


P 3u 

_& E 

pg  3x 


h (T  -T  )4»  s 

-ES  ZR 


u P u P F(u  -u  ) 

_ is  _ _g_g  3£  + & p . 

spg  3x  Px  3x  Pj 

r r . u*\  r r u2  i 

-4Ech“  - - -4E  __Ji+uu 

Pjl  8 2j  p^g  2 g pj 


(A-9) 


The  solid  phase  energy  equation  is  expressed  In  a similar  form  to 


Eq*  CA — 8) | i.e. 9 


<■,»-«%  - - f»„  - » + £ |ppa-*)CEp  + »J/2>J 


p.(l-*)Su_  E_  + 


_R  + -& 


2 

‘ - r Echera  - (a-j 

c p 2 


Note  the  signs  on  the  terms  Fu^S  and  $S  are  opposite  those  in  the  gas 

phase  energy  equation  because  the  overall  mixture  equation  (gas  + solid) 

cannot  have  such  source  terms.  Generally  (as  shown  above),  one  would  have 

to  define  H = E + P /p  , that  is,  the  solid-phase  density,  p and  the 
P P g P r ' P 

gas-phase  pressure,  P is  employed.  Since  the  particles  are  incompressible, 

g 

we  now  drop  the  tern  P /p  in  the  next-to-the-last  term  of  Eq.  (A-10) . 

g P 

To  simplify  Eq.  (A— 10) , one  again  multiplies  the  solid-phase  con- 
tinuity Eq.  (A-3)  by  (Ep  + u*/2)  and  multiplies  the  solid-phase  momentum 
Eq.  (A-7)  by  Up.  Both  of  these  latter  two  equations  are  subtracted  from 
Eq.  (A-10),  and  assuming  q (external  heating)  is  neglected,  one  arrives  at 


- u “ + 
p 3x 


P^  3x 


r . * 
+ — IE  + EChem]  + 

ft  1 n n J ft 


(A-U) 


diem 

Here  Ep  represents  the  latent  heat  of  pyrolsis  for  the  solid. 


CAPTIONS : Figures  1-8 


Figure  1. 

Figure  2. 

Figure  3 
Figure  4 
Figure  5 

Figure  6 

Figure  7 

Figure  8 


Pressure  distribution  during  flamespreading  in  an  initially 
packed  bed  of  small  particles  of  solid  propellant.  (Figs.  2-5 
show  variation  in  other  flow  parameters  for  the  same  case  as 
calculated  for  Fig.  1.) 

Locus  of  ignition  (flame)  front  and  pressure  front,  the  latter 
derived  from  P vs  x distribution. 

Gas  and  particle  velocity  distribution  history. 

Gas  and  particle  temperature  distribution  history. 

Porosity  (gas  volume  fraction)  distribution  history.  The  (*) 
indicates  the  location  of  the  ignition  front. 

Ignition  (flame)  front  locus  for  three  different  initial  spherical 
particle  sizes.  The  insert  shows  at  one  (arbitrary)  time  the 
peak  pressure  in  the  bed  as  a function  of  the  initial  size. 

The  pressure-time  variation  at  six  different  locations  in  the 
bed;  Eq.  (6b)  replaces  Eq.  (6a). 

Locus  of  ignition  (flame)  front  and  the  pressure  front  (derived 
from  observations  of  P(t)  shown  in  Fig.  7)  indicating  rapid 
flame spreading  accerelation  which  may  lead  to  DDT. 


T 


Location  (inches) 

X — - 

Figure  1 Pressure  distribution  history  predicted  for  a three  inch 
long  bed  with  an  initial  solids  loading  of  60  percent. 
[Predictions  for  Figs.  1 - 5 utilized  the  second-form 

of  the  particle  energy  equation.) 


Time  ( fisec ) 

Figure  2 Flame  front  and  pressure  front  loci  for  the  pressure  distributions 
indicated  in  Fig.  1. 


Figure  7 


The  pressure- time  variations  at  six  different  locations. 
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Figure  8 Flame  front  locus  and  pressure  wave  front,  indicating  conditions 
for  possible  DDT. 
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Abs:ra:t 

The  Opper.hein  CLOUD  prccr a.-,  has  been  adapted  to 
study  the  possibility  of  initiating  detonation  in  a 
reactive  mixture  that  contains  a spherically  sym- 
retric  local  region  of  enhanced  reactivity.  The 
induction  tone  for  the  reactive  mixture  is  modeled 
by  an  Arrhenius  kinetic  law.  The  effect  of  the  en- 
hancement is  modeled  by  changing  the  following  fac- 
tors: the  activation  energy  for  the  induction 
tone,  the  spatial  distribution  in  the  enhanced  re- 
activity region,  the  pre-extor.ential  coefficient 
and  the  maximum  allowable  rate  of  energy  release. 
All  these  variables  have  been  found  to  affect  the 
ability  of  the  localized  rapid  reaction  in  the 
central  region  to  trigger  detonation  in  the  sur- 
rounding cloud. 


Introduction 


In  1970  Meyer  et  al . di 
for  initiation  of  detonatior 
taal  shock  heating  is  not  t> 
mechanism.  Recently  Knystau 
shown  experimentally  that  a 
of  enhanced  reactivity  can  c 
of  detonation  without  the  i: 
strong  shock  in  the  system. 


iscovered  a mechanism 
r.  in  which  direct  ini- 
r.e  important  triggering 
utas  et^  al^. 2,3  have  also 
properly  shaped  region 
cause  direct  initiation 
r.itial  presence  of  a 


This  phenomenon  of  direct  initiation  of  de- 
tonation by  localized  enhanced  reactivity  is  very 
important  to  our  understanding  of  the  mechanisms  of 
blast  wave  production  in  unconfined  vapor  cloud 
explosions.  It  is  now  well  known  that  extremely 
high  effective  burning  velocities  must  be  reached 
if  a damaging  blast  wave  is  to  be  generated  by  a 
premixed,  unconfined  cloud  that  is  burning.  Up 
until  the  discovery  of  this  particular  initiation 
phenomenon  the  mechanism  of  transition  to  such  ex- 
tremely high  burning  velocities  or  even  to  detona- 
tion under  conditions  of  partial  confinement  was 
mot  understood.  It  was  obvious  fsom  experimental 
evidence  (the  accidents  themselves)  that  transi- 
tions were  occurring  under  situations  where  the 
confinements  were  not  sufficient  to  first  produce 
strong  shock  waves.  It  currently  appears  that  the 
direct  initiation  by  localized  enhanced  reactivity 
that  Knystautas  and  his  co-workers  have  observed  is 
the  most  viable  mechanism  for  the  production  of 
extremely  high  burning  velocities  or  transition  to 
detonation  in  situations  where  only  partial  con- 
finement is  available. 

It  is  also  well  known  that  constant  volume  ex- 
plosion followed  by  release  of  the  high  pressure 
contents  cannot  generate  shock  waves  that  are 
strong  enough  to  induce  exothermic  chemical  re- 
actions in  hydrocarbon-air  mixtures.4  The  shock 
wave  temperatures  are  just  too  low.  However,  one 
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of  the  mechanisms  proposed  by  Knystautas  et^  aK 
for  the  production  of  detonation  is  a chemical 
reaction-gas  dynamic  process  in  which  the  intense 
chemical  reaction  in  a central  region  generates 
pressure  waves  that  interact  with  a neighboring 
region  having  a spatial  gradient  of  chemical  re- 
activity. Lee  and  co-workers  have  observed  this 
effect  experimentally.  Under  these  conditions  they 
postulated  that  the  waves  enhance  the  reactivity  in 
the  gradient  region  and  a reactive  wave  can  actually 
accelerate  to  detonation  velocity. 

This  paper  presents  the  results  of  a numerical 
study  to  test  the  possibility  that  such  an  accelera- 
tion mechanism  is  indeed  available  in  a completely 
unconfined  situation.  In  the  idealization  used 
here  the  geometry  is  strictly  spherical  and  the 
entire  system  is  assumed  to  be  heated  by  a simple 
thermal  explosion  process  having  hydrocarbon-air 
Arrhenius  kinetics  as  derived  from  shock  tube 
measurements.  In  the  neighborhood  of  the  origin 
the  system  is  arbitrarily  given  enhanced  re- 
activity without  increasing  the  local  heat  of  com- 
bustion. With  these  conditions,  the  CL0UD^»<>  pro- 
gram is  used  to  calculate  the  nonsteady  develop- 
ment of  the  flow  field  for  various  assumed  profiles 
for  the  central  region  of  enhanced  reactivity. 

Mathematical  Formulation 

Our  model  for  this  study  is  heat  addition  to  a 
perfect  gas  with  constant  gamma  throughout.  The 
fluid  is  initially  uniform  end  at  zero  time  the 
region  of  enhanced  chemical  reactivity  is  created 
with  spherical  symmetry  and  without  disturbing  the 
background  gas.  Ke  choose  to  formulate  the  prob- 
lem in  Lagrangian  form,  for  which  the  relevant 
equations  are: 

3v  v 3(ur^) 

3T=-  V 0) 
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where  r is  the  position  coordinate,  y,  v,  u,  e,  p 
and  q,  respectively,  are  the  specific  heat  ratio, 
specific  volume,  radial  velocity,  internal  energy, 
pressure,  and  rate  of  energy  addition,  and  j = 0, 
1,  2 for  planar,  cylindrical  and  spherical 
geometries. 
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A'e  scale  the  problem  by  introducing  rQ,  the 
outer  radius  of  the  region  of  enhanced  reactivity, 

introduce 
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core  over  which  the  enhanced  reactivity  is  uniform. 
In  the  transition  region  where  8 < R i 1,  the 
function  is  given  by 

<t>  (R)  = [cos(3miJi)  - 9cos(tt4i)  ♦ 8]/16  (15) 


i - q = — a-  t . R = r/r  and  t = t/t  . Kith 
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these  quantities,  Eqs.  (l)-(5)  become 
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The  energy  addition  term,  Q,  is  assumed  to  be 
given  by 

Q = — — (11) 

y p v dt  ' 
ro  o 

where  AE  is  the  heat  released  by  a particular 
fuel-oxidizer  system  per  unit  volume  v0,  and  X is 
the  reaction  coordinate  with  range  0 i X < = 1. 
The  ratio  AE/pQV0  is  related  to  the  value  of  M^j 
by  the  expression 

M*  = ! 
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The  reaction  rate  ~ is  assumed  to  be  given  by 
the  Arrhenius  law 


= A SItUE  exp  L I*] 

d-  Ao  n-1  exp  l E J ’ 


where  E*  is  the  dimensionless  activation  energy,  n 
is  the  order  of  the  reaction  and  A is  the  pre- 
exponential factor  which  is  assumed  to  depend  upon 
both  the  reactive  system  and  the  enhanced  chemical 
reactivity  of  the  central  region  as  follows: 


Aq  = Aj  (1  ♦ A) 


0 < R < 8 


= Aj  (1  ♦ A$(R))  8 < R < 1 


1 < R < » 


A is  the  dimensionless  preexponential  factor 
for  the  reactive  system  without  enhancement,  A is 
an  amplifying  factor  which  models  the  effect  of 
enhanced  reactivity,  8(<1)  is  the  radius  of  the 


This  function  blends  very  smoothly  with  the  func- 
tions for  the  adjacent  regions,  for  not  only  do 
their  values  match,  but  their  first,  second  and 
third  derivatives  are  zero  at  the  end  points  and 
therefore  also  match. 

With  the  assumption  of  initial  uniformity,  the 
initial  conditions  are: 


k = 1 ....  .m 


R(0)  = Rfc 
P(Rk,0)  = 1 
V(Rk,0)  = 1 
U(Rk,0)  = 0 
E(Rk,0)  = 2.5 
X(Rk,0)  = 0 


Further,  the  boundary  conditions  at  R = 0 for 
t > 0 are 

UCO.t)  = 0 
|£(0,t)  = 0 

i to.*)  * ° 

§f(0,T)  = 0 

Sr(0'T)  = 0 


The  disturbance  caused  by  the  energy  release 
will  be  propagated  outward  and  the  leading  edge  of 
this  disturbance  will  be  bounded  by  an  acoustic 
wave  or  a shock  wave.  Ahead  of  this  lead  wave,  the 
variables  will  remain  at  their  initial  values  but 
inside  they  change  and  their  values  must  be  com- 
puted. This  is  true  even  though  the  entire  field 
is  assumed  to  be  reactive  with  Arrhenius  kinetics 
(Eq.  13)  because  at  t = 0 we  assume  that  X - 0 
everywhere  and  the  only  cases  of  interest  ^o  this 
study  are  those  when  X «<  1,  unless  the  region  has 
enhanced  reactivity  due  to  some  external  physical 
phenomenon  such  as  flame  folding,  or  the  region  has 
been  shock  heated  to  a condition  that  triggers 
rapid  reaction.  Essentially  this  means  that  the 
nature  of  our  approach  to  the  problem  in  effect 
renders  the  surroundings  nonreactive  during  the 
observation  time.  This  will  be  discussed  further 
in  the  results  section. 

Computation 

The  computer  program  for  the  calculations  is 
based  upon  the  CLOUD  program  of  Oppenheim^  as  sub- 
sequently modified  by  Adamczyk.  The  latter  simply 
added  an  energy  source  term  to  the  energy  equation 
but  retained  the  other  features  of  the  original 


2 


f 


CLOJU  program,  including  the  finite  difference  in- 
tegration scheme  and  artificial  viscosity  to  smear 
out  shock  fronts  or  steep  velocity  gradients.  Of 
course,  we  replaced  Adanczyk's  energy  addition  law 
by  the  Arrhenius  law  detailed  in  Eqs.  (11)  and 
(13). 

Our  first  attempts  at  running  this  program  pro- 
duced very  long  run  times  which  were  caused  by 
extremely  small  time  increments  resulting  from  the 
Courant  stability  requirement  applied  to  the  steep 
gradients  that  developed.  To  shorten  the  run  time 
we  subsequently  introduced  a minimum  time  increment 
which  we  would  allow  and,  as  one  would  expect, 
stability  problems  appeared.  This  instability 
would  manifest  itself  when  new  radii  or  R values 
calculated  by  using  Eq.  (10)  would  not  agree  with 
radii  calculated  from  the  new  volume  or  V values 
based  upon  mass  conservation.  The  problem  would 
a.  ise  when  the  rate  of  energy  addition  had  large 
spatial  gradients,  which  in  turn  caused  the  flow 
velocity,  U,  to  change  rapidly,  so  that  the  average 
value  of  U used  in  the  program  did  not  advance  R 
sufficiently. 

This  instability  problem  was  overcome  by  adding 
one  predictor-corrector  type  iteration  step  to  the 
calculation.  The  order  of  procedure  in  the  CLOUD 
program  is  to  use  Eq.  (7)  and  then  Eq.  (10)  to 
advance  the  velocity  and  position  at  half  points  in 
the  space-time  grid,  then  to  use  Eq.  (6)  to  ad- 
vance the  specific  volume  at  integer  points  in  the 
space-time  grid.  After  advancing  the  specific 
volume  at  point  i,  we  recalculated  the  value  of 
position  at  i ♦ 1/2  by  using  mass  conservation 
with  the  tentative  value  of  specific  volume,  then 
recalculated  the  specific  volume  using  Eq.  (6) 
from  this  corrected  value  of  position. 

Once  the  minimum  time  increment  is  reached  in 
the  calculation,  subsequent  calculations  proceed 
at  this  time  increment  unless  the  time  increment 
as  calculated  from  the  Courant  stability  criterion 
again  becomes  larger.  Ke  have  completed  perhaps 
20  computer  runs  without  difficulty,  and  we  be- 
lieve that  in  some  runs  that  this  minimum  time 
increment  may  have  been  as  much  as  ten  times  the 
value  permitted  by  the  Courant  condition. 

The  Role  of  the  Parameters 

As  formulated,  there  are  eight  dimensionless 
parameters  which  appear  in  the  problem,  namely,  n, 
Y,  6,  !E/p0v0,  E*.  A,  Aj  and  Qq,  where  is  the 
maximum  energy  addition  rate.  The  parameter  n is 
the  order  of  the  reaction  and  has  a characteristic 
value  for  a specific  system  of  reactants,  .he 
specific  heat  ratio,  y,  is  specific  for  each  sys- 
tem. The  parameter  8 defines  the  extent  of  the 
flat  central  core  for  the  region  of  enhanced  re- 
activity and  therefore  the  relative  steepness  of 
the  distribution  in  the  transition  region.  The 
parameter  AE/p0v0  represents  a dimensionless  com- 
bustion energy  and,  as  mentioned  earlier,  de- 
termines the  value  of  shock  Mach  number  which  would 
finally  be  attained  if  a successful  transition  to 
a Chapman- Jouguet  detonation  occurs.  The  dimen- 
sionless activation  energy,  E* , determines  the 
relative  ease  of  reaction  with  temperature  and  is 
linked  to  A and  Aj,  for  the  larger  the  value  of 
E*,  then  the  larger  the  product  of  A and  Aj  must 
be  to  start  the  reaction.  The  amplification  fac- 
tor, A,  essentially  describes  the  "effectiveness" 
of  the  enhanced  reactivity  in  producing  the 


reaction  in  the  core  region. 

The  remaining  two  parameters,  Aj  and  are 
both  essentially  the  ratio  of  two  time  scales,  with 
one  time  scale  always  being  related  to  the  time  for 
a signal  to  propagate  across  the  core  and  the  other 
being  a characteristic  time  for  reaction.  The  role 
of  Aj  is  to  determine  the  incubation  time  in  our 
model  while  ^ determines  whether  transition  to 
detonation  will  occur,  since  another  way  to  inter- 
pret its  meaning  is  that  it  represents  the  ratio  of 
maximum  rate  of  energy  input  to  the  rate  at  which 
energy  can  be  propagated  outward  through  acoustic 
waves.  This  interpretation  implies  a critical 
value  of  Q,,,  below  which  transition  will  not  occur 
and  above  which  it  will  occur.  Further,  if  Aj  is 
too  small,  the  required  critical  may  be  greater 
than  the  system  is  capable  of  producing  and  hence 
implies  that  there  is  a critical  value  of  Aj , say 
AJ,  below  which  the  transition  to  detonation  can- 
not occur.  Consequently,  for  a given  system  with 
localized  enhanced  reactivity  there  is  a minimum 
core  size  which  is  necessary  in  order  to  achieve 
transition  to  detonation.  Further,  for  Aj  values 
larger  than  AJ,  the  system  is  not  actually  con- 
strained as  we  required  in  the  calculations  and 
thus  would  achieve  values  of  larger  than  the 
critical  values  required  for  transition  to  de- 
tonation. These  considerations  are  summarized 
qualitatively  in  Fig.  1. 


Fig.  1 Qualitative  behavior  of  Aj  vs  Qm  for  the 

critical  line  and  maximum  energy  rate  line. 

Numerical  Results 

Our  computational  efforts  have  mainly  been 
based  on  the  following  values  for  the  parameters: 
n=2,  y-1.4,  8=.2,  AE/povo=20,  E*=168  and  A=108.  We 
have  assumed  values  of  Aj  on  the  order  of  1018  and 
computed  the  behavior  for  various  values  of 
which  are  set  by  the  program,  after  we  provide  a 
multiplier  according  to  how  large  we  want  Qjj,  to  be 
compared  to  the  maximum  energy  addition  rate  when 
the  minimum  time  increment  is  determined.  The  E* 
value  is  appropriate  for  propane  or  ethylene  which 
have  activation  energies  on  the  order  of  40  kcal/g 
mole  and  the  AE/p0vQ  value  at  constant  volume  will 
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produce  a pressure  rise  to  initial  pressure  ratio 
equal  to  eight  or  an  Mj-j  value  of  approximately 
5.43. 

In  Fig.  2 the  pressure-position  behavior  at 
various  times  is  plotted  for  = 14.9  and 
Aj  = 6.96  x 10  . Clearly  this  case  failed,  for 

the  lead  wave  reached  a pressure  rise  ratio  of 
approximately  1.6  near  the  edge  of  the  enhance- 
ment zone  at  R = 1 and  then  started  to  decay.  The 
numerical  output  reveals  that  outside  the  enhanced 
region  the  reaction  takes  place  so  slowly  that  the 
energy  addition  does  not  overcome  the  effect  of  the 
spherical  divergence,  so  the  wave  decays. 
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Fig.  2 Pressure  rise  ratio  vs  dimensionless  radius 
at  different  dimensionless  times  for  a ^ 
value  at  which  transition  to  detonation 
failed.  Top,  during  time  interval  while 
enhanced  region  reacts;  Bottom,  during  time 
interval  while  very  little  reaction  occurs 
anywhere.  Aj=6.96  x 1018.  (^=14.9. 

This  result  should  be  contrasted  with  Fig.  3 
which  shows  the  pressure-position  behavior  at 
various  times  for  = 28.9  and  the  same  Aj.  Now, 
while  there  appears  to  be  a momentary  hesitancy  as 
the  lead  wave  passes  R = 1,  the  strength  continues 
to  grow,  albeit  rather  slowly,  for  at  the  end  time 
shown,  the  pressure  rise  ratio  is  only  about  251  of 
the  pressure  rise  ratio  for  a Chapman-Jouguet  wave 
for  this  system.  The  important  fact  is  that  the 
gas  is  reacting  outside  the  region  of  enhancement 
to  increase  the  strength  of  the  lead  wave. 
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Fig.  3 Pressure  rise  ratio  vs  dimensionless  radius 
at  different  dimensionless  times  for  a 
value  at  which  transition  to  detonation  is 
occurring  very  slowly.  Top,  during  time 
interval  while  enhanced  region  reacts; 
Bottom,  during  time  interval  while  reaction 
occurs  outside  enhanced  region. 

Aj=6.96  x 1018.  (fcn-28.9. 

We  have  similar  results  also  for  Aj=6. ObxlO1^. 
For  §,,,=26. 8 the  wave  failed,  while  for  <^=53 . 6 the 
wave  strength  was  still  growing  at  the  end  of  the 
computer  run. 

Additional  features  we  noticed  in  the  computer 
outputs  include  the  following:  At  the  larger  value 
of  Aj  the  uniform  core  reacts  essentially  at  con- 
stant volume,  but  at  the  smaller  value  of  A,  the 
reaction  moves  out  progressively  in  the  uniform 
core.  For  both  values  of  At,  when  is  low  the 
reaction  region  is  very  diffuse  and  may  involve 
ma:ny  "cells",  perhaps  as  many  as  one  hundred  while 
for  larger  ^ the  reaction  region  is  much  thinner, 
being  of  the  order  of  ten  cells.  Also  the  larger 
the  value  of  Q,,,,  the  more  rapidly  the  strength  of 
the  lead  wave  approaches  the  value  for  a CJ  wave. 
This  behavior  is  shown,  for  example,  by  these  data 
for  cases  with  Aj=6.96xl0*8.  With  Qm=28.9  the 
maximum  pressure  ratio  has  risen  to  only  9.72  at 
t=1 .89  while  for  Qm=149.3  the  maximum  pressure 
ratio  has  risen  to  19.1  at  t=1.07. 

For  Aj=6.96xl018,  the  "incubation"  time,  de- 
fined as  the  time  for  a cell  to  reach  a pressure 
ratio  of  1.08,  is  .275  while  for  Aj=6.96xl07  the 
incubation  time  is  4.42,  more  than  15  times 
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We  close  by  noting  that  typical  values  of  X in 
the  region  outside  the  leading  wave  are  of  the 
enter  of  10'*^  at  the  end  of  the  computer  run. 

These  values  are  so  stall  that  the  surrounding  re- 
gi:n  effectively  is  r.onreactive  during  the  time 
period  of  interest. 

Discussion 

Obviously,  the  results  presented  here  are 
largely  preliminary.  Nonetheless,  these  results 
confirm  that  direct  initiation  with  localized  en- 
hanced reactivity  is  possible  in  our  model  for 
certain  combinations  of  parameter  values  and  a 
particular  distribution  in  the  enhancement  transi- 
tion region.  Further,  they  seer,  to  confirm  that 
the  transition  region  is  necessary  but  indicate 
nothing  concerning  the  optimum  shape  of  that  re- 
gion. Clearly,  they  show  that  some  minimum  wave 
strength  is  needed  at  the  edge  of  the  enhancement 
region  if  the  reaction  wave  is  to  continue  to  grow 
ir.  amplitude  as  it  moves  outside  the  enhancement 
region.  Thus,  even  at  this  early  stage  of  de- 
velopment, this  concept  shows  that  limiting  the 
maximum  local  rate  of  energy  addition  is  a possible 
method  to  prevent  transition  to  detonation  by 
1 realized  enhanced  reactivity. 

Based  on  these  preliminary  results,  it  is  con- 
templated that  a complete  parametric  study  will  be 
performed. 
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